#include "mecafluide.hpp"

// variational formulation language
//#include <life/lifevf/vf.hpp>

//#include <mecafluide_cl.cpp>

namespace Feel
{
namespace FSI
{
    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J, vector_ptrtype& R)
    {

        using namespace Feel::vf;

        if (M_LookComment) std::cout <<"[NavierStokes] : updateJacobian start \n";
        boost::timer btime;
        btime.restart();

        this->changeRepository();

        //mesh_ptrtype mesh = Xh->mesh();
        auto mesh = this->mesh();
        auto Xh = this->functionSpace();


        element_fluid_type U( Xh, "u" ); U = *X;
        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = U.template element<0>();
        element_fluid_0_type v = V.template element<0>();
        element_fluid_1_type p = U.template element<1>();
        element_fluid_1_type q = V.template element<1>();

        //--------------------------------------------------------------------------------------------------//


        //Tenseur des deformations
        auto deft = 0.5*(gradt(u)+trans(gradt(u)) );
        auto def = 0.5*(grad(v)+trans(grad(v)) );
        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmat = -idt(p)*Id + 2*M_mu*deft;
        // Tenseur des contraintes (test)
        auto Sigma = -id(q)*Id + 2*M_mu*def;

        auto clDef = MECAFLUIDE_CL(M_time);
        //--------------------------------------------------------------------------------------------------//
        //--------------------------------------------------------------------------------------------------//
        //--------------------------------------------------------------------------------------------------//

        form2( Xh,Xh, J , _init=true );

        //terme non lineaire
        bool M_isNLTerme=true;
        if (M_isNLTerme)
            {
                form2( Xh,Xh, J ) +=
                    integrate (elements(mesh), _Q<3*Order>(),
                               //M_rho*trans(id(v))*(gradv(u))*idt(u));
                               M_rho*trans(idt(u))*(gradv(u))*id(v));
                form2( Xh,Xh, J )  +=
                    integrate (elements(mesh), _Q<3*Order>(),
                               //M_rho*trans(id(v))*(gradt(u)*idv(u)) );
                               M_rho*trans(idv(u))*(gradt(u)*id(v)) );

                if (M_isMoveDomain)
                    {
                        form2( Xh,Xh, J )  +=
                            integrate (elements(mesh), _Q<3*Order>(),
                                       -M_rho*trans(idv(*M_mesh_velocity))*(gradt(u))*id(v) );
                    }
            }

#if 0
        // sigma : grad(v) sur Omega
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh), _Q<2*(Order)>(),
                       + trace(Sigmat*trans(grad(v))) );
#else
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh),
                       + 2*M_mu*trace(trans(deft)*grad(v)) );
                       //+ 2*M_mu*trace(deft*trans(grad(v))) );
                       //+ 2*M_mu*trace(trans(deft)*def) );
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh),
                       - idt(p)*div(v) );
#endif


#if 1
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh),// _Q<2*(Order)>(),
                       - M_rho*divv(*M_mesh_velocity)*trans(idt(u))*id(v) );
#endif


        // terme assurant l'incompressibilité
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh), _Q<2*(Order-1)>(),
                       + divt(u)*id(q) );


        //--------------------------------------------------------------------------------------------------//



        if (M_weakCL)
            {
                ForEachCL( clDef, cl::dirichlet_vec,
                           form2( Xh, Xh, J ) +=
                           integrate( markedfaces(mesh,PhysicalName), //_Q<2*Order>(),
                                      //- trans(Sigmat*N())*id(v)
                                      - 2*M_mu*trans(deft*N())*id(v) + trans(idt(p)*N())*id(v)
                                      + penalbc*trans(idt(u))*id(v)/hFace() ) );

                ForEachCL( clDef, cl::paroi_mobile,
                           form2( Xh, Xh, J ) +=
                           integrate( markedfaces(mesh,PhysicalName), //_Q<2*Order>(),
                                      - trans(Sigmat*N())*id(v)
                                      //- trans(Sigma*N())*idt(u) // sym
                                      //- 2*M_mu*trans(deft*N())*id(v) + trans(idt(p)*N())*id(v)
                                      + penalbc*trans(idt(u))*id(v)/hFace() ) );

            }
        // terme necessaire pour condition forte avec petsc
        form2( Xh, Xh, J ) +=
            integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
                       +0*div(v)*idt(p)+0*divt(u)*id(q) +0*idt(p)*id(q) );

        //bool M_isStat=true;
        if (!M_isStat)
            {
                form2( Xh, Xh, J ) +=
                    integrate( elements(mesh), _Q< 2*Order >(),
                               + M_rho*trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );
            }



        double penalisationBeta=0.05;
        auto beta_N = val(trans(N())*M_rho*(idv(u)-idv(*M_mesh_velocity)));
        double order_scaling = math::pow( double(Order), 3.5 );
        form2( Xh, Xh, J )
            += integrate( internalfaces(Xh->mesh()), //typename MyIm< 4*uOrder-2 >::type(),
                          val(penalisationBeta + vf::pow(hFace(),2.0)*( beta_N*beta_N ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )/order_scaling
                          );

#if 0
        int stabilisation = this->application()->vm()["stab"].template as<int>();
        if (stabilisation==1) {
            //terme de stabilisation
            form2( Xh, Xh, J ) +=
                integrate( internalfaces(mesh), _Q< 3*Order >(),
                           + 0.025*vf::pow(hFace(),2.0)*( trans(idt(u))*N() )*(trans(jump(grad(u)))*jump(grad(v)) )  );
            //+ 0.025*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
            //            + (0.05+vf::pow(hFace(),2.0))*( (trans(idv(u))*N())*( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
            //( (trans(idt(u))*N())*( trans(id(u))*N() )) );

            form2( Xh, Xh, J ) +=
                integrate( elements(mesh), _Q< 3*Order >(),
                           + 0.000*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
        }
#endif


        J->close();

        if (!M_weakCL)
            {
                auto dirElem = this->getDirichletElement();
                ForEachCL( clDef,cl::paroi_mobile,
                           form2( Xh, Xh, J ) +=
                           on( markedfaces(Xh->mesh(), Xh->mesh()->markerName(PhysicalName)) , *dirElem, R, vf::vec(vf::cst(0.),vf::cst(0.)) ) );

                ForEachCL( clDef,cl::dirichlet_vec,
                           form2( Xh, Xh, J ) +=
                           on( markedfaces(Xh->mesh(), Xh->mesh()->markerName(PhysicalName)) , *dirElem, R, vf::vec(vf::cst(0.),vf::cst(0.)) ) );

            }

        double timeElaps = btime.elapsed();
        if (M_LookComment) std::cout <<"[NavierStokes] : updateJacobian finish in " << timeElaps << "s\n";

    }

} // end namespace NavierStokes

} // Feel


// instantiation
template class Feel::FSI::NavierStokes<MECAFLUIDE_DIM,MECAFLUIDE_ORDER_VELOCITY,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;

